Model of light collimation by photonic crystal surface modes 
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We propose a quantitative model explaining the mechanism of light collimation by leaky surface 
modes that propagate on a corrugated surface around the output of a photonic crystal waveguide. 
The dispersion relation of these modes is determined for a number of surface terminations. Analytical 
results obtained on the basis of the model are compared to those of rigorous numerical simulations. 
Maximum collimation is shown to occur at frequency values corresponding to excitation of surface 
modes whose wave number retains a nonzero real part. 
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I. INTRODUCTION 

One of the problems hindering wider commercial ap- 
plication of photonic crystals (PCs) is the difficulty in 
coupling PC waveguides to conventional dielectric wave- 
guides or optical fibers. A possible solution consists in 
tapering the waveguide so as to achieve better coupling 
with the fiber; this has been the subject of a number of 
publications, e.g., Refs.[l||2|. Recently, Moreno et al. 3 and 
Kramper et al. 4 independently suggested that collimation 
of the light emitted by a waveguide (beaming) could also 
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occur due to excitation of surface modes in the proximity 
of the waveguide exit. Based on the earlier discovery of 
a similar effect in metallic structures supporting surface 
plasmons^ the idea has been expanded in several articles 
following the original papers j&LSaS n 1Q i 1 1 

Moreno et al. propose a simple qualitative theory to ex- 
plain the novel effect 4 The radiation reaching the wave- 
guide outlet can couple to surface modes; if the surface 
around the outlet is corrugated (i.e., modulated with pe- 
riod different from that of the underlying crystal), its 
eigenstates become 'leaky', since energy is emitted as the 
radiation scatters at the perturbed surface cells. Un- 
der appropriate conditions, the scattered waves interfere 
constructively along the surface normal, thus producing 
a collimated beam. According to Ref. y, this constructive 
interference takes place for surface modes of wave vector 
k x = (shifted to the first Brillouin zone of the surface), 
with the phase difference between two successive scatter- 
ers equal to an integer multiple of 2tt. 

The purpose of this work is to formulate a quantita- 
tive model of the beaming effect in PCs, taking explicitly 
into account the imaginary component of the leaky mode 
wave vector. The model predictions are tested against 
results of numerical simulations. We also show that, for 
practically realizable PCs, maximum beaming occurs for 
surface modes with Re k x 7^ 0. 



II. MODEL 

The system to be considered is depicted in Fig. [TJ Ex- 
cited by a source at the waveguide input (left), a guided 
mode propagates towards the crystal surface (right). On 
reaching the output, the radiation is partially reflected, 
partially emitted directly into free space, and the remain- 
der excites leaky modes propagating upwards and down- 
wards along the surface corrugated with period A. Either 
side of the waveguide comprises N surface unit cells. Our 
aim is to calculate the radiation intensity @(0) defined as 



FIG. 1: A waveguide embedded in a PC with corrugated sur- 
face. 



0(0) = lim rS r (r,0), 
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where S r (r, 0) denotes the radial component of the time- 
averaged Poynting vector at the point specified by the 
polar coordinates (r, 6). 

Without losing generality, we restrict our attention to 
E polarization (with the electric field parallel to the cylin- 
der axes). Since all sources of the electromagnetic field 
are located in the z < halfspace, the crystal can be 
regarded as an aperture antenna, with radiation pattern 
proportional to the Fourier transform of the electric field 
distribution at the z = axis: 12 

E(r,0) -ye^-^ 4 ) -^/(/c sin<9)cos6>, (2a) 
ykor 

where 

f(k sin 6) = -= / E y (x,0)e- ikoXBine dx, (2b) 

ko = uj/c denotes the free-space wave number and r is 
assumed to be large compared to the system dimensions. 
Since the crystal is symmetric with respect to the z axis, 



E y (x,0) = E y (— x,0) and 

f(k smO) =f+(k sm6)+f+(-k sm6), (3) 
where / + (fcosin#) is defined as 

/+ (k sin 6) = —= / E y (x,0)e- ikoxsine dx, (4) 
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The field E y (x,0) consists basically of three major com- 
ponents: the beam stemming directly from the waveguide 
outlet [Ey g (x, 0)], the leaky mode propagating along the 
corrugated surface [i^ urf (x, 0)], and the residual radia- 
tion extending past the crystal boundaries [E y es (x, 0)]. 
In this section, we focus on the surface mode contribu- 
tion, assuming E y (x,0) = £^ urf (:r, 0); the effect of the 
other components is discussed in Section HVl 

Neglecting the fringe effects at the surface boundaries, 
we can apply the Bloch theorem to the field related to 
the leaky mode. This yields 



| x | _ |) e iMM-d/2) if o < |a;| - \ < NA, 







otherwise. 
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Function u(x) is yl-periodic and k x = k' x + \k" (&" > 0) denotes the leaky mode wave number. We assume the 
waveguide output has effective width <i, which is the lattice constant of the underlying PC, and do not consider this 
area to belong to the crystal surface. With Eq. ([5]) substituted into Eq. (|4|), we get through integration: 
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where the structure factor F(kosmQ) is defined as 
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i(x)e [ ^- kosine)x dx. (7) 



Being a periodic function, u(x) can be Fourier-expanded: 

oo 

u(x) = J2 w«e 27rinx/yl , (8) 
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resulting in the following form of the formula for 
F(k sm6): 

F(k Sin 0) = —= V U n / e Kk x +2nn/A-k sin 9)x dx 
V 27r n J <> 
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In the first approximation, which is often used in analyt- 
ical treatment of leaky- wave antennas^ only the term 
with denominator of the smallest magnitude (the zeroth 
term in the case of near-zero k x ) needs to be kept in the 
above sum. Substitution of this approximate structure 
factor into Eq. (|6|) leads to 
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From the Maxwell equations it can be shown that in vac- 
uum 
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Sr(r,0) = --\E(r,e)\ 2 with Z = J (11) 
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(9) Therefore, using Eqs. (j2a|). dHU)— (JTTD, and the defi- 
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nition (pQ), we arrive at 
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The unknown coefficient uo can be determined from 
the principle of energy conservation. Consider a semi- 
infinite section of the crystal surface along which the 
leaky mode in question propagates. By a procedure anal- 
ogous to that followed above, we obtain the radiation in- 
tensity generated by the leaky wave in this configuration: 
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Integrated over the interval [— f , f ], this yields the total 
power radiated into free space, which must be equal to 
the power Pq exciting the leaky mode. Therefore, 
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The integral in the above equation can be evaluated an- 
alytically. Thus, 
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where 
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By substituting Eq. ([T5|) into Eq. ([12]), we obtain the 
final formula for the radiation intensity in the system 
shown in Fig. [TJ normalized to the accepted power P 09 
and expressed solely in terms of the leaky surface mode 
parameters: 
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III. NUMERICAL DETERMINATION 
OF SURFACE MODES 

To apply the above-discussed model to a specific pho- 
tonic surface, e.g., for the determination of the frequency 
most conducive to beaming, it is necessary to calculate 
the dispersion relation of the modes supported by the sur- 
face. In this section we shall briefly outline the method 
we employed for this purpose. 

In our approach, we consider a semi-infinite PC with 
possible surface reconstruction. The whole system is 
divided into three parts: the homogeneous region, the 
surface, and the underlying semi-infinite, but otherwise 
ideal, photonic crystal. The electromagnetic field in the 
homogeneous material is represented as a Rayleigh ex- 
pansion, i.e., a linear combination of discrete plane waves, 



whereas in the semi-infinite crystal the field is expanded 
into the eigenmodes of the corresponding infinite struc- 
ture (a procedure suggested by Istrate et ali&). Since 
we are searching for states leaking energy away from 
the surface, in both regions we only consider waves that 
propagate or decay in this direction; the precise rules 
of choosing these waves are specified in the Appendix. 
The complex band structure necessary to find the field 
representation in the PC is calculated by the differential 
method (see Ref. [l5| for details). 

The fields in the homogeneous region and in the PC 
are linked by the scattering matrix— of the surface layer, 
which provides the necessary boundary conditions. This 
leads to a homogeneous system of linear equations, which 
must have a non- trivial solution for surface states to exist. 
The search for surface modes is thus reduced to a search 
for roots of the determinant of a matrix dependent on 
k' x , k" and uj. 

IV. RESULTS 

A. Surface mode dispersion 

In the following we shall focus on crystals of the type 
shown in Fig. [H considering a truncated square lattice 
of dielectric cylinders of permittivity e = 11.56 embed- 
ded in vacuum, with bulk cylinder radius 0.18<i, surface 
cylinder radius 0.09<i, surface corrugation period A = 2d, 
and three values of corrugation depth z corr : 0, — O.ld 
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FIG. 2: Crystal A surface mode dispersion curve. The 
hatched areas represent bulk bands, and the grey triangle de- 
notes the bound-wave region, in which all spatial harmonics 
are evanescent in vacuum for k" — 0. 



and — 0.3d (the minus sign indicating that the perturbed 
cylinders are shifted towards the crystal). For future ref- 
erence, we denote these three crystals with letters A, £>, 
and C, respectively. Fig. [2] presents the dispersion curve 
of surface modes supported by crystal A (in the extended 
Brillouin zone scheme), while in Fig. [3] the dispersion 
curves corresponding to crystals B and C are plotted 
in the vicinity of the center of their first Brillouin zone. 
All these dispersion relations have been determined by 
the method described in the previous section. Only the 
modes with k x > are shown. 

In the frequency range uo < 0.409 x 2itc/d crystal A 
supports proper (non- leaky, i.e., with k x = 0) surface 
waves, since this fragment of the dispersion curve lies in 
the bounded- wave region (shaded in Fig. [2]), in which all 
spatial harmonics are evanescent in vacuum. Figure (H[a) 
shows the field map of a sample mode from this part of 
the dispersion curve. The range 0.409 < uod/2itc < 0.438 
corresponds to a stop band, where the surface mode 
wave number takes values n/d + \k x \ see Fig. \MJ®) for 
the field map of a typical stop-band state. Character- 
istic for periodic structures, the occurrence of a stop 
band at the Brillouin zone boundary has been observed 
in surface-mode dispersion relations of periodic dielec- 
tric waveguides embedded in homogeneous media j 16 i 17 In 
contrast, the shape of the dispersion curve at k' x > n/d 
stems directly from the presence of the underlying PC, 
i.e., from the periodicity of the 'substrate'. Although 
the dispersion curve remains within the bounded-wave 
region, the surface mode wave number retains a large 
imaginary part, since the wave leaks energy into the crys- 
tal, as can be observed in Fig. ffl(c), showing the field 
magnitude of the mode labeled III in Fig. 03 This field 
structure is analogous to that of leaky modes propagating 
in periodic waveguides adjacent to homogeneous media: 
the wave amplitude in the PC grows as z — > — oo. As 
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FIG. 3: Crystal B and C surface mode dispersion curves 
(dashed and dotted lines, respectively) in the vicinity of the 
center of the first Brillouin zone. The solid line represents a 
segment of the crystal A dispersion curve from Fig. [2] shifted 
by Ak x = —ir/d and shown for comparison. 



observed in early studies of leaky waves, 18 this behavior 
is not unphysical, since in real systems the field extent is 
limited by the location of the source exciting the leaky 
wave. 

Let us proceed to the case of nonzero corrugation 
(crystals B and C). Doubling the surface period re- 
sults in the first Brillouin zone folding away to the range 
[— 7r/2d, 7r/2d]; consequently, at frequency values above 
0.25 x 2irc/d at least one spatial harmonic is radiative 
and the originally bound surface states become leaky. 
As shown in Fig. O the surface mode dispersion curve 
smoothes out, shifting towards negative k' x and larger 
positive k x , as the corrugation depth increases. Inter- 
estingly, at frequency values above 0.438 x 2itc/d the 
substrate PC supports bulk states characterized by es- 
sentially imaginary k x (k x ~ ik x > 0.66i x 2k/ d) and 
purely real k z . In consequence, in all three crystals con- 
sidered here, the leaky modes from the immediate vicin- 
ity of the k' x = line are propagative in the — z direc- 
tion. Therefore, this line serves as a boundary between 
leaky modes exponentially decaying (k' x < 0) and grow- 
ing (k x > 0) inside the substrate. Figures Sfd)-(f) show 
maps of the field corresponding to points IV- VI on the 
dispersion curve of crystal C. 



B. Light collimation: frequency dependence 

Let us now proceed to the analysis of beaming itself. 
Figure [5] presents the electric field magnitude calculated 
by the multiple-scattering (MS) method (see Ref. Ojjl for 
details), with geometry parameters and frequency value 
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FIG. 4: (a)-(f) Maps of electric field magnitude for surface 
modes labeled I-VI in Fig. [3] The frequency and wave number 
values are expressed in units of 2nc/d and 27r/d, respectively. 



conducive to directional emission. In Fig. El the frequency 
dependence of the radiation intensity ^> s i m (# = 0) calcu- 
lated by the MS method is compared to that obtained on 
the basis of our model, for crystals B and C with N = 9 
and N = 15 corrugations. In these MS simulations we 
consider the system depicted in Fig. [TJ with the wave- 
guide n wg = 12 cylinders long and cladding n c i a d = 5 
cylinders wide; a waveguide mode is excited by a point 
source near the inlet. The results depicted in Fig. [6] 
clearly show that our model reproduces the basic feature 
of the effect in question, i.e., the existence of a distinct 
transmission maximum at a well-defined frequency value. 
The relative height of the @ S im(6 = 0) curve peaks for dif- 
ferent crystals is rendered reasonably well too. There are 
some visible differencies between the model predictions 
and the simulation results, though. Most notably, the 
absolute maxima in the theoretical curves are shifted by 
approximately 0.002 x 2irc/d to the right with respect 
to those found numerically. Furthermore, at frequency 
values smaller (larger) than those corresponding to the 
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FIG. 6: Frequency dependence of the collimated beam inten- 
sity for crystal B (a) and C (b). The left and right y axes 
refer to the results of MS calculations (lines) and the model 
predictions (circles), respectively. 



main peaks, the model generally predicts radiation inten- 
sity values several times smaller (larger) than those calcu- 
lated numerically. Possible causes of these discrepancies 
are analyzed in the following subsection by scrutinizing 
the angular dependence of the radiation intensity. 



C. Light collimation: angular dependence 

Figure [71(a) shows the angular dependence of the ra- 
diation intensity ^ s i m (#) calculated by the MS method 
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FIG. 7: (a) Angular dependence of the radiation inten- 
sity for crystal C with N = 9 corrugations at frequency 
u — 0.400 x 2nc/d, calculated by the MS method. Solid black 
line: total radiation intensity; dotted line: contribution of the 
waveguide outlet region, calculated by Eq. J2} with the inte- 
gration interval restricted to \x\ < d/2; dashed line: contribu- 
tion of the surface region, with d/2 < \x\ < NA + d/2; grey 
line: contribution of the crystal exterior, with \x\ > NA-\-d/2. 
(b) Angular dependence of the radiation intensity for the same 
crystal and frequency value, calculated on the basis of our 
model with the surface field expansion [Eq. J8}] truncated to 
a single harmonic (solid line) and to 17 harmonics (dashed 
and dotted lines). The data plotted with the dotted line 
result from computations taking into account surface wave 
reflections at the crystal boundaries. 



for crystal C with N = 9 corrugations at frequency 
ou = 0.400 x 27Tc/d. To help evaluate the relative impor- 
tance of the three regions mentioned in Section [II] — the 
waveguide outlet, the corrugated surface and the crystal 
exterior — we have also plotted separately the contribu- 
tions of sections \x\ < d/2, d/2 < \x\ < NA + d/2 and 
\x\ > NA + d/2, calculated by the aperture formula ([2]) 
with the values of E y (x,0) obtained by the MS method. 
The radiation intensity produced by each of these regions 
will be labelled an ^ ^s?m> respectively. 

The graph shows that while surface modes play the 
most important part and are responsible for the forma- 
tion of the major lobes, the radiation coming directly 
from the waveguide outlet has an impact as well. In 
particular, it makes the main beams shift by several de- 
grees towards the surface normal. As a result, they begin 
to overlap at frequency values below those predicted by 
leaky- wave considerations alone. This explains the dis- 



FIG. 8: Same as Fig. but at frequency oo = 0.420 x 2nc/d. 
Since at this frequency virtually no power reaches the crystal 
boundaries, ^sfm(^) an d the reflected surface wave contribu- 
tion are not calculated. 

crepancy in position of the absolute maxima between the 
theoretical and the numerical curves in Fig. [6], mentioned 
in the previous subsection. 

The radiation pattern in the large-angle (0 > 15°) re- 
gion is also affected by the field stemming from the outlet: 
the destructive interference of this field with that emitted 
from the corrugated surface causes an offset of the side- 
lobe positions. The contribution of the residual radia- 
tion, however, remains negligible throughout the angular 
range covered by the plots in Fig. 

Let us compare the results presented in Fig. 0(a) with 
the predictions of our model, plotted in Fig. 0(b). The 
solid curve has been calculated with the approximate 
leaky mode structure factor [i.e., with the series in Eq. J5J) 
truncated to a single term, as discussed in Section [II] ; 
the dashed line results from calculations with the full 
structure factor (with amplitudes of the individual har- 
monics computed numerically by the method outlined in 
Section Hill) . In both curves, the main peaks occur at 
#o ~ 11°, in good agreement with the ^irn curve plotted 
in Fig. 0(a). However, the field structure at angles far 
from 0o is seen to depend strongly on the structure fac- 
tor; in general, the radiation intensity obtained with the 
one-harmonic approximation is much smaller than that 
calculated without this simplification. This may be the 
reason why at low frequency values the model-predicted 
<P(0 = 0) value is very small compared to that resulting 
from MS calculations, since for these frequency values 
the surface normal lies far from the direction = 0q. 
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Another difference between the theoretical plots in 
Fig. 0(b) and the curve is a distinct dip of the latter 
at 9 « 6°. This turns out to follow from surface wave re- 
flections occurring at the crystal boundaries, which have 
been neglected in our model, but are of some importance 
for crystals with shallow corrugations (involving weakly 
leaky modes) or of small size; see Fig. [9] The dotted 
curve in Fig. 0(b) shows the radiation intensity angular 
dependence after taking these reflections into account, 
with the reflection coefficient at the surface termination 
calculated by the method outlined in Ref. [2(J. Evidently, 
this curve is in excellent agreement with the results of 
MS simulations. 

Lastly, let us consider the frequency uo = 0.420 x 27rc/d, 
which lies to the right of the main peak in the ^ S im(# = 0) 
curve shown in Fig. [61(b), and for which the model- 
predicted radiation intensity at 6 = is too large. 
Fig. [81(a), a counterpart of Fig. 0(a), shows that at 
this frequency value the radiation stemming from the 
waveguide output has magnitude larger than in the low- 
frequency case analyzed previously. However, since in 
the small-angle region it is out-of-phase with the field 
produced by the leaky modes, the total radiation inten- 
sity becomes significantly smaller than it would be with- 
out the direct beam from the waveguide outlet. In addi- 
tion, as indicated by Fig. [HI neglecting higher harmonics 
in the surface field expansion used in the model-based 
calculations leads to some overestimation of the leaky- 
mode-induced radiation intensity at = 0. Together, 
these two factors explain the theory-versus-simulation 
discrepancy observed in Fig. [6] in the frequency range 
u > 0.41 x 2nc/d. 



V. DISCUSSION 

As indicated in Fig. [HI maximum beam collimation oc- 
curs at frequency value of around 0.41 x 2nc/d, which, 
according to Fig. El corresponds to surface modes with 
k' x « -0.01 x 27r/d (crystal B) and k' x « -0.025 x 2n/d 
(crystal C) rather than to modes from the center of the 
Brillouin zone (k x = 0). Incidentally, the k' x = modes, 
being delocalized (propagative inside the crystal), could 
not be responsible for beaming. However, the nonzero 
real part of the wave vector of the surface modes for 
which maximum beaming is observed is easily explained 
on the basis of the model discussed in Section [TTJ. It is 
a consequence of the competition between the tendency 
to reduce \k' x \ in order to obtain better phase-matching 
of waves emitted from individual unit cells, and, on the 
other hand, the negative effect of too large a k x on the 
effective length of the 'grating'. Since in the negative- 
k' x region a decrease in \k' x \ is always accompanied by 
an increase in k x: the most intensive beaming occurs for 
moderate ('balanced') values of both parameters. 

The model also sheds new light on the fact that in 
negatively corrugated crystals the frequency value corre- 
sponding to maximum beaming lies remarkably close to 



that of the unperturbed surface mode from the Brillouin 
zone center. This proves to be a resultant of two opposing 
effects. It has been pointed out 10 that bringing surface 
cylinders closer to the bulk crystal causes a blueshift of 
the surface mode dispersion relation, due to a decrease of 
the fraction of electromagnetic energy contained within 
the dielectric. However, since beaming occurs at nonzero 
k' x values, the shift starts from an initial frequency value 
lower than that of the Brillouin-zone-center mode. As a 
consequence, the resultant optimum beaming frequency 
is close to the original frequency of the k x = mode. 

The accuracy of the model could be considerably im- 
proved by taking into account the radiation emitted di- 
rectly from the waveguide outlet. This, however, would 
require a detailed investigation of the interactions be- 
tween the outlet and the surface cylinders in its immedi- 
ate vicinity, as these interactions determine the amount 
of power transferred to surface modes and that emit- 
ted directly into free space. An analytical formulation 
of these effects seems hardly feasible, though. 



VI. CONCLUSIONS 

We have presented a quantitative analysis of a model 
explaining the effect of surface corrugation on the colli- 
mation of radiation leaving the outlet of a photonic crys- 
tal waveguide, on the basis of the dispersion relation of 
leaky modes supported by the corrugated surface. The 
dispersion relation has been calculated and discussed for 
a number of surface terminations. The model has been 
shown to reproduce the basic features of the investigated 
effect, and the significance of the factors left out of ac- 
count has been evaluated. Besides clarifying the con- 
ditions necessary for optimum beaming, the model also 
explains the relative insensitivity of the maximum colli- 
mation frequency value to the degree of surface modula- 
tion. We believe our results will contribute to a deeper 
understanding of the physical grounds of the beaming 
effect. 
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APPENDIX: THE CHOICE OF BASIS STATES 

In this section we specify the selection rules for the 
states to be used as the expansion basis for fields in the 
homogeneous and periodically modulated system regions 
considered in Section UTTl Let us begin with the homoge- 
neous region. Here, the electric field can be written as a 
Rayleigh expansion, i.e., a linear combination of spatial 



FIG. 9: Electric field magnitude in the surface region of crystal B with N = 9 corrugations at frequency uo — 0.405 x 27rc/d. 
An interference pattern resulting from surface wave reflections at the crystal boundaries is clearly visible. 



harmonics: 

E(x,z) = ^A n e^ kxnX+kznZ \ (A.l) 

n 

where k xn = k x + 2irn/A, and the condition k zn — k\ — 
k 2 n holds; k x and ko are fixed. For real k Xl the obvious 
choice for the sign of k zn is 

k zn > if k 2 zn > 0, (A.2a) 
k zn /i>0 if *4<0, (A.2b) 

so that the propagating harmonics (|A.2a|) carry energy in 
the +z direction, and the evanescent ones (|A.2b|) decay 
with z — > oo. When k x is complex, the sign of k zn should 
be chosen so as to assure analytical continuity with the 
k x G M case, i.e., 

Re /c 2n > if Re ^ n > 0, (A.3a) 
Im k zn > if Re k 2 zn < 0. (A.3b) 

It is easy to show that these rules can be combined into 

Rek zn +Imk zn > 0, (A.4) 



in accordance with Refs. fl6ll2ll . 

In the semi-infinite PC, the field can be expanded in 
the crystal eigenstates corresponding to the fixed k x and 
fco- For real k x , these eigenstates are either purely prop- 
agative in the z direction, with k z G R, or purely evanes- 
cent, with k z = ik z or k z = ±7r/d+ifc", where k z G R. In 
the former case, the states to be included in the expan- 
sion are those with negative z component of their energy 
flux vector £, and in the latter case, those decaying when 
z approaches — oo, i.e., those with k" z < 0. When k x is 
allowed to be complex, the k z component of a crystal 
eigenstate wave vector can take arbitrary complex values 
too. However, the eigenmodes can still be classified as 
'essentially propagative', with k z = k z + ik z fulfilling the 
condition 

\K\<\k' z \ and \K\<7r/d-\k' z \, (A.5) 

and 'essentially evanescent' otherwise. States of these 
two families to be included in the expansion should then 
be selected according to the £ • z < and k z < rules, 
respectively. 
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